Notas previas:

  • Todos los archivos se encuentran en /home/alumno14/analisis_datos_omicos/tarea/
  • El código ejecutado en QIIME2 se encuentra en metagenomic/metagenomica_codigo_QIIME.sh
  • Los archivos generados tras el análisis metagenómicos se metagenomic/analisis.
  • El archivo del análisis de R de metagenómica se encuentra en metagenomic/Graficas.Rmd
  • El manifiesto y los metadatos de metagenómica se encuentran en metagenomic/manifest.csv y metagenomic/metadata.tsv
  • El código ejecutado para el alineamiento del RNA-seq se encuentra en rnaseq/rnaseq.sh
  • El pipeline de GALAXY se encuentra en rnaseq/Galaxy.html
  • El análisis en profundidad del RNA-seq se encuentra en rnaseq/‘RNASeq ADO.Rmd’ y rnaseq/RNASeq-ADO.html
  • El análisis en profundidad del scRNA-seq se encuentra en scrnaseq/‘Single Cell.Rmd’ y scrnaseq/Single-Cell.html

Introducción y justificación

Los objetivos específicos de la práctica se dividen en varias fases analíticas:

  • El análisis metagenómico explora cómo la microbiota inicial puede influir en la respuesta del sistema nervioso central a la inflamación sistémica, buscando correlaciones entre la microbiota y las respuestas a la sepsis.

  • El análisis de RNAseq busca identificar mecanismos moleculares y celulares que contribuyen a la neurodegeneración y neuroinflamación en un modelo murino de sepsis, con un enfoque en las diferencias de género.

  • El análisis de scRNAseq profundizará en la respuesta celular específica a la sepsis, centrando la atención en células como los macrófagos microgliales y células linfoides para entender las dinámicas de población celular.

  • Finalmente, en futuros estudios se planea una integración multiómica de los datos de RNAseq, scRNAseq, y metagenómica para descubrir interacciones complejas y posibles dianas terapéuticas, utilizando datos proporcionados por el grupo de Oftalmología Experimental de la Universidad de Murcia.

Metagenómica

Se va a realizar un análisis comparativo del metagenoma de ratones criados en la Universidad de Murcia frente a aquellos criados en la Universidad de Leuven, Bélgica. Este estudio se motiva por observaciones preliminares que indican que los ratones de Murcia exhiben una mayor sensibilidad y una respuesta inflamatoria exacerbada al tratamiento con LPS (lipopolisacárido), a diferencia de los ratones de Leuven, que muestran una reacción atenuada o no lo presentan.

La secuenciación metagenómica de estas muestras se llevó a cabo utilizando la tecnología Ion Torrent. Para el análisis y visualización de los datos obtenidos, se empleará QIIME 2 (Quantitative Insights Into Microbial Ecology 2), una herramienta de software de código abierto especialmente diseñada para estudiar la ecología microbiana a partir de secuencias de ADN. Este enfoque permite identificar diferencias en la composición de la microbiota que podrían explicar las variaciones en la respuesta al LPS entre las cohortes murinas de ambas universidades.

Metodología

Importación y control de calidad

El primer paso crucial en el uso de QIIME 2 es la importación adecuada de los archivos FASTQ. Para nuestro estudio, disponemos de 19 archivos FASTQ correspondientes a 19 muestras diferentes, que ya están demultiplexadas, según se detalla en el archivo metadata.tsv. Para la importación, utilizamos el formato Casava 1.8, y optamos por la estrategia de importación mediante un archivo de manifiesto. Este archivo de texto simple es fundamental para que QIIME 2 pueda interpretar correctamente la estructura de los datos de entrada. En él se mapean las muestras a sus respectivos archivos de secuencias, especificando la ubicación exacta de los archivos y la manera en que deben ser procesados. El manifiesto debe incluir, al menos, dos columnas: una para el identificador de cada muestra y otra para las rutas a los archivos de secuencias (consultar manifest.csv).

Una consideración esencial al importar archivos Casava 1.8 mediante manifiesto en QIIME2 es que los nombres de los archivos no deben contener el carácter “_”. Este carácter puede causar problemas durante el procesamiento posterior, ya que QIIME 2 podría no reconocer el archivo correctamente.

El manifiesto se genera a partir de las siguiente linea de comandos:

echo "sample-id,absolute-filepath,direction" > manifest.csv
ls *.fastq | awk -v pwd="$PWD/" '{sub(/\.fastq$/, "", $1); print $1 "," pwd $1 ".fastq,forward"}' >> manifest.csv

Lo siguiente es ejecutar la línea de comandos para importar los ficheros FASTQ a qiime para el procesamiento de los mismos. Para visualizar la calidad de las lecturas de los ficheros FASTQ se ejecuta seguidamente la siguiente linea:

qiime tools import \
  --type 'SampleData[SequencesWithQuality]' \
  --input-path manifest.csv \
  --output-path demux.qza \
  --input-format SingleEndFastqManifestPhred33
  
qiime demux summarize \
  --i-data demux.qza \
  --o-visualization demux.qzv

La Figura 1 de la izquierda muestra la calidad promedio de las lecturas de nuestras muestras. Observamos que la calidad se mantiene en un promedio de 30 para las primeras 280 bases, pero luego se vuelve inestable. Por esta razón, decidimos excluir esta sección inestable de las secuencias.

Dado que nuestra herramienta está configurada para procesar lecturas de secuenciación de Illumina, necesitamos realizar algunos ajustes en el protocolo de procesamiento. Conforme a Catozzi et al. 2019, es recomendable recortar las primeras 15-20 bases de cada secuencia para eliminar los primers utilizados en la secuenciación. Además, seguimos la recomendación de Salipante et al., 2014, de descartar lecturas con menos de 100 bases.

Para realizar estos ajustes, empleamos un workflow en Galaxy. Iniciamos con la herramienta “FASTQ Trimmer by column” (Galaxy Version 1.1.5) para eliminar las 20 primeras bases y recortar las secuencias hasta la base 280. Posteriormente, utilizamos “Filter FASTQ reads by quality score and length” (Galaxy Version 1.1.5) para filtrar aquellas secuencias que superen las 100 bases de longitud. Finalmente, estas secuencias filtradas se importan de nuevo a QIIME2 para su análisis posterior, como se ilustra en la Figura 1B.

Figura 1: Control de calidad de los FASTQ. En la imagen de la izquierda podemos ver el resultado antes del procesado. A la derecha con un recorte de las 20 primeras pb y de las últimas, además del filtrado de lecturas mayores de 100 pb.
Figura 1: Control de calidad de los FASTQ. En la imagen de la izquierda podemos ver el resultado antes del procesado. A la derecha con un recorte de las 20 primeras pb y de las últimas, además del filtrado de lecturas mayores de 100 pb.

Inferir features

Existen dos estrategias para calcular la diversidad de una muestra. Mediante Unidades Taxonómicas Operativas (OTUs) o mediante inferencia de Variantes de Secuencia de Amplicones (ASVs). Las OTUs se agrupan por similitud de secuencia, generalmente con un umbral del 97%, lo que significa que las secuencias que son al menos un 97% similares se consideran la misma OTU. Las ASVs no se agrupan por similitud; cada variante única de secuencia se trata como una entidad distinta.

Según el manual de deblur de QIIME2, que genera OTUs, este comando está diseñado específicamente para datos generados a partir de protocolos de amplicones 16S en plataformas Illumina: “This mode of operation should only be used when data were generated from a 16S amplicon protocol on an Illumina platform”. Siguiendo las recomendaciones de Catozzi et al. y diversas fuentes en blogs, se aconseja emplear la herramienta DADA2 para este proceso. A diferencia del método tradicional que genera OTUs, DADA2 permite la creación de ASVs, que proporcionan una resolución taxonómica más precisa. Para implementar esto, ejecutamos lo siguiente:

qiime dada2 denoise-single \
 --i-demultiplexed-seqs demux.qza \
 --p-trunc-len 0 \
 --p-n-threads 0 \
 --o-representative-sequences rep-seqs-dada.qza \
 --o-table table-dada.qza \
 --o-denoising-stats denoising-stats.qza
Tabla 1. Tabla que muestra la longitud y el número de secuencias encontradas en nuestras muestras
Statistic count min max mean range std
Value 6117 25 260 230.386 235 34.2403

Vemos que se han inferido un total de 6117 secuencias distintas con un tamaño medio de 230 pb (Tabla 1).

Para elegir un sampling-depth adecuado, primero necesitamos examinar la distribución de las profundidades de secuenciación de nuestras muestras. Generalmente hay que elegir una profundidad que mantenga el máximo número de muestras posibles mientras se sacrifica lo menos posible en términos de número de secuencias. La mediana de la frecuencia por muestra que encontramos que es de 101.820 features. Sin embargo, decidimos escoger una sampling depth de 30.000 features, recogiendo así el 25,25% de las features en el 78,95% de las muestras.

Arbol filogenético

La generación de un arbol filogenético es fundamental para entender las relaciones evolutivas entre las secuencias de especies en nuestras muestras. La generación de un árbol filogenético se basa en diferentes métricas como la diversidad filogenética de Faith, UniFrac ponderado y no ponderado. El método se basa en un alineamiento múltiple de las secuencias para generar la filogenia en base a la similitud de las mismas. Este paso es esencial para el posterior análisis de métricas de diversidad alfa y beta.

Ya tenemos todo lo necesario para hacer el árbol filogenético. Esto se genera mediante los siguientes comandos:

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree-dada.qza \
  --i-table table-dada.qza \
  --p-sampling-depth 30000 \
  --m-metadata-file metadata.tsv \
  --output-dir core-metrics-results \
  --p-n-jobs-or-threads auto \
  --parallel

Para observar si hemos recogido suficiente diversidad en nuestras muestras realizamos un análisis de Alpha rarefacción (Figura 2). Para ello ejecutamos lo siguiente:

qiime diversity alpha-rarefaction \
  --i-table table-dada.qza \
  --i-phylogeny rooted-tree-dada.qza \
  --p-max-depth 30000 \
  --m-metadata-file metadata.tsv \
  --o-visualization alpha-rarefaction-dada.qzv
**Figura 2:** *Gráfico de rarefación. Índice shannon de cada muestra según la profundidad de secuenciación.*

Figura 2: Gráfico de rarefación. Índice shannon de cada muestra según la profundidad de secuenciación.

El análisis de la gráfica de rarefacción alfa (Figura 2) nos indica que todas las muestras alcanzan un máximo en el índice de Shannon, lo cual sugiere que incrementar el número de secuencias probablemente no revelará un número significativo de especies nuevas. Esto indica que por más profundidad de secuenciación no encontraremos especies nuevas. Además, se observan variaciones en la diversidad alfa entre las diferentes muestras; algunas registran valores más elevados en el índice de Shannon que otras. Estas diferencias podrían atribuirse tanto a variaciones reales en la biodiversidad de las comunidades muestreadas como a la eficacia del muestreo o de los métodos de secuenciación empleados.

Análisis

En el análisis de metagenómica, es habitual calcular diversos índices de diversidad que proporcionan una medida de la abundancia y riqueza de especies en cada muestra. Estos índices facilitan la comparación de estas características entre diferentes muestras o grupos de muestras.

La diversidad alfa se refiere a la riqueza de especies dentro de una unidad de muestreo específica, que puede ser una muestra individual o un grupo de muestras bajo la misma condición. Para medirla, se utilizan varias métricas como el índice de Shannon, la riqueza de especies observadas, equidad de Simpson, diversidad inversa de Simpson y la diversidad filogenética de Faith, entre otras, que consideran tanto la presencia y ausencia de especies como su abundancia relativa.

Por otro lado, la diversidad beta evalúa la variación o cambio de diversidad entre diferentes ambientes o condiciones, utilizando métricas de distancia que reflejan la disimilitud entre muestras, siendo las más comunes la distancia de Bray-Curtis, el Índice de Jaccard y Unifrac, para determinar cuán similares o diferentes son las muestras entre sí.

Para obtener todas las métricas para estudiar en profundidad ejecutamos la siguiente línea de comandos:

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree-dada.qza \
  --i-table table-dada.qza \
  --p-sampling-depth 30000 \
  --m-metadata-file metadata.tsv \
  --output-dir core-metrics-results \
  --p-n-jobs-or-threads auto \
  --parallel

Diversidad Alfa

Con el código anterior se generan los ficheros para observar estas métricas y poder hacer estudios estadísticos al respecto. Como vemos en la Figura 3 y la Figura 4, en la comparación entre grupos de sexo o entre universidades, respectivamente, no observamos ninguna diferencia en la diversidad o abundancia recogidas en las muestras. Si nos fijamos en profundidad, cada universidad tiene unos outliers (2 en la KU y 3 en la UMU) con una diversidad menor que el resto. Esto podría ser debido a un fallo en la secuenciación o que estas muestras tienen algo distinto al resto de muestras.

**Figura 3:** *Medidas de estudio de diversidad alfa en las muestras distinguiendo por Sexo. De izquierda a derecha: Indice Shannon, Abundancia, Faith e indice de equidad.*

Figura 3: Medidas de estudio de diversidad alfa en las muestras distinguiendo por Sexo. De izquierda a derecha: Indice Shannon, Abundancia, Faith e indice de equidad.

**Figura 4:** *Medidas de estudio de diversidad alfa en las muestras distinguiendo por Universidad. De izquierda a derecha: Indice Shannon, Abundancia, Faith e indice de equidad.*

Figura 4: Medidas de estudio de diversidad alfa en las muestras distinguiendo por Universidad. De izquierda a derecha: Indice Shannon, Abundancia, Faith e indice de equidad.

Diversidad beta

Vamos ahora a estudiar la diversidad beta para ver si se parecen las especies de microorganismos entre grupos. Tanto la distancia de Bray-Curtis (Figura 5) como la distancia de Jaccard (Figura 6) son medidas para medir la disimilitud entre muestras. Podemos observar que en el PCoA de ambas distancias se diferencian claramente 4 poblaciones dentro de nuestras muestras y estan claramente definidas por las Universidades, ya que el sexo no es un factor diferencial para separar las poblaciones .

Figura 5: PCA plot de la distancia Bray-Curtis para estipular la beta diversidad entre las muestras. Las muestras de color morado pertenecen a la Universidad de Leuven y las de amarillo a la UMU. Los puntos circulares son muestras de sexo femenino y los cuadrados, masculinos.

Figura 6: PCA plot de la distancia Jaccard para estipular la beta diversidad entre las muestras. Las muestras de color morado pertenecen a la Universidad de Leuven y las de amarillo a la UMU. Los puntos circulares son muestras de sexo femenino y los cuadrados, masculinos.

Si hacemos un estudio estadístico de los datos, observamos que con ambas métricas (Bray-Curtis y Jaccard) hay una diferencia significativa entre las muestras de cada Universidad (Tabla 2 y Tabla 3).

Tabla 2: Resultados del Pairwise permanova con respecto a la distancia Bray-Curtis entre universidades
Group.1 Group.2 Sample.size Permutations pseudo.F p.value q.value Metric
KU UMU 15 999 3.565576 0.001 0.001 Bray-Curtis
Tabla 3: Resultados del Pairwise permanova con respecto a la distancia Jaccard entre universidades
Group.1 Group.2 Sample.size Permutations pseudo.F p.value q.value Metric
KU UMU 15 999 2.801261 0.001 0.001 Jaccard

Composición taxonómica

Cada ASV identificada previamente necesita ser asignada a un taxón específico, y para esta tarea, en QIIME2, se puede emplear un clasificador Bayesiano que ha sido preentrenado. Dicho clasificador se ha derivado de la base de datos Greengenes de OTUs al 99%, la cual ha sido recortada para ajustarse a lecturas de 250 bases correspondientes a la región V4 del ARNr 16S.

Para ello ejecutamos el siguiente código:

qiime feature-classifier classify-sklearn \
  --i-classifier gg-2022-10-backbone.v4.nb.qza \ 
  --i-reads rep-seqs-dada.qza \
  --o-classification taxonomy-dada.qza \
  --p-n-jobs 5 \
  --verbose
  
qiime metadata tabulate \
  --m-input-file taxonomy-dada.qza \
  --o-visualization taxonomy-dada.qzv

qiime taxa barplot \
  --i-table table-dada.qza \
  --i-taxonomy taxonomy-dada.qza \
  --m-metadata-file metadata.tsv \
  --o-visualization taxa-bar-plots-dada.qzv

En este gráfico podemos ver la variedad de especies encontradas en las muestras (Figura 7). Como vemos en la figura de la izquierda, cada muestra tiene una variedad distinta. Cuando comparamos la variedad relativa de las muestras por universidades y sexo vemos que existen poblaciones distintas según el grupo.

**Figura 7:** *Abundancia relativa de taxones por muestra, universidad y sexo.*

Figura 7: Abundancia relativa de taxones por muestra, universidad y sexo.

Análisis diferencial

Para observar si hay diferencias en la composición microbiana entre ratones de las distintas universidades, aplicaremos DESeq2 a los conteos taxonómicos. Esta herramienta realiza ajustes para la profundidad de secuenciación y variabilidad biológica, y mediante modelos estadísticos, nos ayudará a identificar de manera fiable y precisa los taxones cuya abundancia difiere significativamente entre las muestras.

dds <- DESeqDataSetFromMatrix(countData=counts.taxa, 
                              colData=metadata, 
                              design=~University, tidy = TRUE)

dds <- DESeq(dds)
res <- results(dds)
res <- res[order(res$padj),]
rownames(res) <- sub(".*\\.(.*)$", "\\1", rownames(res))
res <- data.frame(res)

Observamos que las especies más significativamente distintas entre universidades se muestran en la tabla siendo la especie Prevotella la más significativa con un log2FoldChange de 10 veces más en la Universidad de Leuven que en la Universidad de Murcia (Tabla 4).

Tabla 4: Tabla con las medidas del análisis diferencial por medio DESeq2 de las 5 especies más diferencialmente distintas entre universidades
baseMean log2FoldChange lfcSE stat pvalue padj
Prevotella 115.710659 10.035202 1.075937 9.326945 1.089666e-20 1.155046e-18
Parasutterella 58.977861 -25.773993 2.961139 -8.704080 3.201618e-18 1.696857e-16
distasonis 60.277236 9.090762 1.081676 8.404334 4.302956e-17 1.520378e-15
massiliensis 5.395228 21.154224 2.972465 7.116728 1.105195e-12 2.928766e-11
caecigallinarius 129.895009 10.186081 1.581145 6.442218 1.177398e-10 2.496085e-09
veroni666751 1459.347035 -9.119939 1.439587 -6.335106 2.371776e-10 4.190137e-09

En el volcano plot podemos ver las especies más representadas en una Universidad o en otra (Figura 8). Esas mismas especies podrían ser las responsables que los ratones de la Universidad de Leuven se vean menos afectados al tratamiento por LPS y no entren en un estado de sepsis severo.

**Figura 8:** *Volcanoplot que muestra la relación entre el cambio (log2 Fold Change) y la significancia estadística (-log10 p-valor) de especies entre universidades. Las especies más presentes en la Universidad de Leuven están representados en amarillo, los de la Universidad de Murcia en morado y los no significativos en gris. *

Figura 8: Volcanoplot que muestra la relación entre el cambio (log2 Fold Change) y la significancia estadística (-log10 p-valor) de especies entre universidades. Las especies más presentes en la Universidad de Leuven están representados en amarillo, los de la Universidad de Murcia en morado y los no significativos en gris.

RNA-seq

Utilizaremos un modelo murino de sepsis para investigar los mecanismos celulares y moleculares que subyacen a la neurodegeneración y la neuroinflamación. Nos centraremos en componentes del inflamasoma y emplearemos RNAseq para ampliar nuestra comprensión de la respuesta del SNC a la inflamación sistémica y sus diferencias entre géneros, según se describe en estudios previos (Rodríguez-Ramírez et al., 2024).

Metodología

Preparación y Control de Calidad de las Secuencias

Inicialmente, se recibieron doce archivos en formato FASTQ, correspondientes a doce muestras distintas: seis de machos (M) y seis de hembras (F), subdivididas en controles y tratados con LPS.

Para asegurar la calidad de las secuencias, se utilizó la herramienta FASTQC, disponible en el servidor Galaxy del máster de bioinformática. FASTQC evalúa la calidad de las secuencias de datos de secuenciación de alto rendimiento, proporcionando información sobre la calidad del base calling, la distribución de la calidad de las secuencias, el contenido de bases y la presencia de secuencias contaminantes, entre otros (Figura 9A y B).

Figura 9: Análisis de calidad pre-procesamiento mediante FASTQC de la muestra 1MTreated. (A) Calidad de la secuencia por base. (B) Contenido de la secuencia por base.
Figura 9: Análisis de calidad pre-procesamiento mediante FASTQC de la muestra 1MTreated. (A) Calidad de la secuencia por base. (B) Contenido de la secuencia por base.

Aunque la mayoría de las muestras mostraron una calidad aceptable, las 10 primeras y últimas bases presentaron problemas de calidad en varias muestras. Por lo tanto, se decidió recortar estas regiones utilizando Fastq Trimmer, una herramienta que permite recortar secuencias de FASTQ en posiciones específicas. Posteriormente, se realizó una segunda evaluación de calidad con FASTQC para confirmar la mejora en la calidad de las secuencias (Figura 10 A y B).

Figura 10: Evaluación de calidad de la muestra de secuenciación denominada 1MTreated tras el procesamiento con FastQ Trimmer. (A) Calidad de la secuencia por base. (B) Contenido de la secuencia por base.
Figura 10: Evaluación de calidad de la muestra de secuenciación denominada 1MTreated tras el procesamiento con FastQ Trimmer. (A) Calidad de la secuencia por base. (B) Contenido de la secuencia por base.

Alineamiento de las Secuencias

Para el alineamiento de las secuencias, se utilizó HISAT2, una herramienta de alineamiento que permite el mapeo eficiente de secuencias de ARN contra un genoma de referencia. En este caso, se empleó el genoma de ratón (versión mm9.fa.gz), disponible en el servidor Galaxy. HISAT2 es particularmente útil para el alineamiento de RNAseq por su capacidad de identificar empalmes de exones y adaptarse a la variabilidad genética (Figura suplementaria 1).

Evaluación de la Calidad del Alineamiento

Posteriormente, se evaluó la calidad del alineamiento mediante Samtools stats, una herramienta que proporciona estadísticas detalladas del alineamiento, incluyendo la proporción de lecturas correctamente alineadas, la cobertura del genoma, y la distribución de las lecturas respecto a regiones genómicas específicas (Figura suplementaria 2).

Cuantificación de la Expresión Génica

Para la cuantificación de la expresión génica, se empleó StringTie, que utiliza los archivos BAM generados por HISAT2 para ensamblar y cuantificar los transcritos. StringTie es capaz de reconstruir con precisión los transcritos y estima las abundancias de estos utilizando un modelo estadístico que permite la comparación entre muestras. Para este proceso, se utilizó el archivo de anotaciones gencode.vM10.annotation.gtf como referencia (Tabla 5).

Tabla 5. Encabezado de la tabla que proporciona una matriz detallada de counts de RNA que enumera los genes (filas) y las muestras analizadas (columnas).
Tabla 5. Encabezado de la tabla que proporciona una matriz detallada de counts de RNA que enumera los genes (filas) y las muestras analizadas (columnas).

Preparación de Datos para Análisis en R

Una vez obtenidos los conteos de expresión de cada muestra mediante StringTie, se descargaron todos los archivos de conteo. Utilicé una función personalizada en Python para concatenar estos archivos en una sola tabla de datos (ver fichero adjunto Merge_Count_tables.py), que luego se utilizó para realizar el análisis estadístico y bioinformático en R.

Carga y Preparación de Datos

El análisis en R comenzó con la importación de los datos de expresión génica combinados usando la función read.csv para cargar los datos desde un archivo CSV. Se realizó una limpieza inicial de datos para asegurar que no hubieran duplicados, utilizando la función aggregate para sumar los conteos en caso de encontrar genes duplicados. Esta preparación de datos aseguraba que cada gen estuviera representado de forma única en el análisis subsiguiente.

Análisis Exploratorio de Datos

Para una primera comprensión de los datos, se emplearon técnicas de análisis exploratorio como el Análisis de Componentes Principales (PCA), que facilita la visualización de las diferencias globales entre las muestras. Utilizando las funciones DESeq para la normalización y plotPCA de la librería DESeq2, se analizaron las principales variaciones en los datos, permitiendo identificar agrupaciones naturales o outliers.

Análisis Diferencial de Expresión Génica

El núcleo del análisis consistió en el uso de DESeq2, un paquete de Bioconductor para análisis de expresión diferencial que utiliza modelos estadísticos basados en la distribución negativa binomial para testar diferencias en la expresión génica. Este análisis se realizó comparando las muestras de tratados versus controles, así como entre machos y hembras, permitiendo identificar genes que mostraban cambios significativos de expresión bajo diferentes condiciones.

Visualización de Resultados

Para la visualización de los resultados del análisis diferencial, se utilizaron varios enfoques:

  • Gráficos de Volcán: Empleando ggplot2, estos gráficos muestran la magnitud de los cambios en la expresión (log2 fold change) versus la significancia estadística (-log10 p-value), facilitando la identificación rápida de genes significativamente regulados.
  • Heatmaps: Utilizando pheatmap, se visualizaron los patrones de expresión de genes significativos, ofreciendo una representación intuitiva de cómo los genes están agrupados en diferentes muestras y condiciones.
  • Diagramas de Venn: Con el paquete ggvenn, se compararon los genes diferencialmente expresados entre diferentes grupos para identificar genes exclusivos y compartidos entre condiciones y sexos.

Enriquecimiento de Genes

Para interpretar los resultados en un contexto biológico más amplio, se realizaron análisis de enriquecimiento utilizando clusterProfiler y org.Mm.eg.db para asociar los genes diferencialmente expresados con rutas biológicas conocidas. Esto incluyó el uso de la función enrichGO para el enriquecimiento de términos de ontología génica y enrichKEGG para identificar rutas metabólicas clave implicadas en las respuestas observadas. Los resultados de estos análisis se visualizaron utilizando enrichplot y pathview para proporcionar interpretaciones gráficas de las rutas y procesos biológicos afectados.

Análisis

Comparación machos controles contra machos enfermas

La comparación inicial se realizó entre los machos controles y los machos tratados con LPS, con el objetivo de identificar diferencias en la expresión génica que podrían elucidar las respuestas del SNC a la inflamación sistémica.

En el análisis de PCA, observamos una anomalía con la muestra de un ratón tratado (2MTreated), que mostró una proximidad notable a las muestras de control (Figura 11A). Esto sugiere una posible falta de respuesta al tratamiento o un error técnico durante la administración del mismo. Al examinar la expresión de genes que mostraron una diferenciación clara en otros ratones tratados, confirmamos que la muestra 2MTreated no seguía el mismo patrón, destacando en particular en el gen Meg1 (maternal expressed gene 3) (Figura 11B). Por lo tanto, se decidió excluir esta muestra del análisis posterior (Figura 11C).

Figura 11: Análisis comparativo de la expresión génica en machos control versus machos tratados con LPS. (A) Representación de análisis de componentes principales (PCA) de datos de RNAseq. (B) CountPlot que refleja que el perfil de expresión del gen Meg3. (C) Representación de análisis de componentes principales (PCA) de datos de RNAseq con la exclusión de la muestra 2MTreated.
Figura 11: Análisis comparativo de la expresión génica en machos control versus machos tratados con LPS. (A) Representación de análisis de componentes principales (PCA) de datos de RNAseq. (B) CountPlot que refleja que el perfil de expresión del gen Meg3. (C) Representación de análisis de componentes principales (PCA) de datos de RNAseq con la exclusión de la muestra 2MTreated.

Una vez ajustado el conjunto de datos, realizamos representaciones gráficas como el volcano plot y el heatmap (Figura 12A y 12B). Estas visualizaciones confirmaron la presencia de múltiples genes diferencialmente expresados entre los machos tratados y los controles, con un patrón de expresión homogéneo en las muestras restantes. Los genes más destacados incluyen Cebpd (CCAAT/enhancer binding protein delta), Lrriq4 (leucine-rich repeats and IQ motif containing 4), Tmem256 (transmembrane protein 256), Ghdc (GH3 domain containing), Meg3 y un pseudogén relacionado con Meg3, todos con cambios significativos de expresión (Figura 12C).

Figura 12: Representaciones gráficas para análisis diferencial de expresión génica. (A) Volcanoplot que muestra la relación entre el cambio en la expresión (log2 Fold Change) y la significancia estadística (-log10 p-valor) de genes entre machos tratados versus control. Los genes significativamente sobreexpresados están representados en amarillo, los subexpresados en morado y los no significativos en gris. (B) Heatmap que ilustra el patrón de expresión de genes en muestras individuales, con los colores oscilando de púrpura (baja expresión) a amarillo (alta expresión). (C) CountPlot de los seis genes más relevantes identificados, Cebpd, Lrriq4, Tmem256, Ghdc, Meg3 y 9630010A21Rik, donde la altura de las barras representa el conteo normalizado de la expresión, diferenciando entre machos control (azul) y tratados (amarillo), resaltando las variaciones sustanciales en la expresión génica inducidas por el tratamiento con LPS.
Figura 12: Representaciones gráficas para análisis diferencial de expresión génica. (A) Volcanoplot que muestra la relación entre el cambio en la expresión (log2 Fold Change) y la significancia estadística (-log10 p-valor) de genes entre machos tratados versus control. Los genes significativamente sobreexpresados están representados en amarillo, los subexpresados en morado y los no significativos en gris. (B) Heatmap que ilustra el patrón de expresión de genes en muestras individuales, con los colores oscilando de púrpura (baja expresión) a amarillo (alta expresión). (C) CountPlot de los seis genes más relevantes identificados, Cebpd, Lrriq4, Tmem256, Ghdc, Meg3 y 9630010A21Rik, donde la altura de las barras representa el conteo normalizado de la expresión, diferenciando entre machos control (azul) y tratados (amarillo), resaltando las variaciones sustanciales en la expresión génica inducidas por el tratamiento con LPS.

En cuanto a las vías biológicas afectadas, el análisis de Gene Ontology señaló la producción de citoquinas por leucocitos mieloides, la producción y regulación de interleucina 1, y la regulación negativa del crecimiento celular o la apoptosis como procesos significativamente alterados. Adicionalmente, el análisis de vías de KEGG identificó la Enfermedad Hepática Alcohólica, la vía de señalización de ErbB, la formación de trampas extracelulares de neutrófilos y la señalización de HIF-1 como rutas destacadas (Figura 13A, B y C).

Figura 13: Análisis de enriquecimiento de vías biológicas en machos tratados frente a control. (A) GSEA Plot que indica enriquecimiento de genes en la producción de citoquinas por leucocitos mieloides, con una curva de enriquecimiento que muestra la distribución de la puntuación de enriquecimiento acumulativo a través del conjunto de genes clasificados. (B) Dotplot de GO:BP (Gene Ontology: Biological Processes) que destaca los procesos biológicos alterados, donde el tamaño del punto indica el recuento de genes y la escala de color representa el p valor, señalando procesos como la producción de interleucina-1 y la regulación de la apoptosis. (C) Dotplot de KEGG (Kyoto Encyclopedia of Genes and Genomes) que muestra las vías afectadas, con el tamaño del punto reflejando la cantidad de genes involucrados y la coloración indicativa del p valor, resaltando vías como la de la enfermedad hepática alcohólica y la señalización de HIF-1.
Figura 13: Análisis de enriquecimiento de vías biológicas en machos tratados frente a control. (A) GSEA Plot que indica enriquecimiento de genes en la producción de citoquinas por leucocitos mieloides, con una curva de enriquecimiento que muestra la distribución de la puntuación de enriquecimiento acumulativo a través del conjunto de genes clasificados. (B) Dotplot de GO:BP (Gene Ontology: Biological Processes) que destaca los procesos biológicos alterados, donde el tamaño del punto indica el recuento de genes y la escala de color representa el p valor, señalando procesos como la producción de interleucina-1 y la regulación de la apoptosis. (C) Dotplot de KEGG (Kyoto Encyclopedia of Genes and Genomes) que muestra las vías afectadas, con el tamaño del punto reflejando la cantidad de genes involucrados y la coloración indicativa del p valor, resaltando vías como la de la enfermedad hepática alcohólica y la señalización de HIF-1.

Estos resultados sugieren que la respuesta a LPS en ratones machos puede estar mediada por alteraciones significativas en la expresión de genes clave y la activación de vías moleculares críticas asociadas con la respuesta inflamatoria y la regulación celular, proporcionando así una base para futuras investigaciones sobre los mecanismos subyacentes de la neuroinflamación y neurodegeneración en contextos de sepsis.

Comparación hembras controles contra hembras enfermas

En nuestro estudio sobre la respuesta inflamatoria inducida por LPS en ratones, la comparación entre las hembras control y las hembras tratadas con LPS reveló hallazgos significativos, destacando diferencias claras en la expresión génica. Utilizando el análisis de PCA, confirmamos la validez de todas las muestras incluidas en el estudio, sin la necesidad de excluir ninguna por comportamiento atípico (Figura 14A).

Figura 14. PCA resultante de la comparación de las muestras hembras control y hembras tratadas. En color azul se representan las hembras tratadas, y en color rojo, las hembras control.
Figura 14. PCA resultante de la comparación de las muestras hembras control y hembras tratadas. En color azul se representan las hembras tratadas, y en color rojo, las hembras control.

El análisis mediante volcano plot demostró un número superior de genes diferencialmente expresados en las hembras en comparación con los machos, acompañado de mayores valores de Fold Change y p-values más bajos, indicando cambios más pronunciados y estadísticamente más significativos (Figura 15A). El heatmap reflejó un comportamiento consistente entre las muestras de ratones control y tratados, con patrones de expresión uniformes para los genes representados (Figura 15B).

Figura 15: Visualizaciones gráficas del análisis de expresión génica en hembras tratadas frente a hembras control. (A) Volcanoplot que muestra diferencias en la expresión génica, con genes sobreexpresados indicados en amarillo, subexpresados en morado y aquellos no significativos en gris, reflejando la magnitud del cambio de expresión (log2 Fold Change) y la significancia estadística (-log10 p-valor). (B) Heatmap que presenta los patrones de expresión génica en las muestras individuales, diferenciando las muestras de control y tratamiento para hembras a través de una gradiente de color que va desde el púrpura (menor expresión) al amarillo (mayor expresión).
Figura 15: Visualizaciones gráficas del análisis de expresión génica en hembras tratadas frente a hembras control. (A) Volcanoplot que muestra diferencias en la expresión génica, con genes sobreexpresados indicados en amarillo, subexpresados en morado y aquellos no significativos en gris, reflejando la magnitud del cambio de expresión (log2 Fold Change) y la significancia estadística (-log10 p-valor). (B) Heatmap que presenta los patrones de expresión génica en las muestras individuales, diferenciando las muestras de control y tratamiento para hembras a través de una gradiente de color que va desde el púrpura (menor expresión) al amarillo (mayor expresión).

Dentro de los seis genes más diferencialmente expresados, encontramos coincidencias con los machos en genes como Meg3, Lrriq4 y Ghdc. Sin embargo, en las hembras también se identificaron genes no anotados como Gm43473 y Gm35339, junto con Mcm10 (minichromosome maintenance 10 replication initiation factor), que son de particular interés debido a sus potenciales roles en la replicación del ADN y la regulación celular (Figura 16).

Figura 16: Análisis de conteo de expresión para los seis genes más diferencialmente expresados en hembras tratadas en comparación con hembras control: Gm42473, Ghdc, Lrriq4, Gm35339, Mcm10, Meg3 y Gm35339. Las gráficas individuales representan cada gen, con puntos morados indicando hembras control y puntos amarillos para hembras enfermas.
Figura 16: Análisis de conteo de expresión para los seis genes más diferencialmente expresados en hembras tratadas en comparación con hembras control: Gm42473, Ghdc, Lrriq4, Gm35339, Mcm10, Meg3 y Gm35339. Las gráficas individuales representan cada gen, con puntos morados indicando hembras control y puntos amarillos para hembras enfermas.

Respecto a las vías biológicas alteradas, observamos la misma actividad en las rutas identificadas en los machos, con la adición significativa de la ruta de regulación de la contracción del músculo liso en las hembras (Figura 17A y B). El análisis de las vías de KEGG resaltó la participación de vías implicadas en varios tipos de cáncer, como el cáncer de páncreas y la leucemia mieloide, así como procesos metabólicos específicos como el metabolismo de la colina y la ruta de señalización AGE-RAGE (Figura 17C).

Figura 17: Análisis de enriquecimiento de vías biológicas basado en datos de expresión génica de hembras tratadas versus control. (A) Dotplot de Gene Ontology Biological Processes (GO:BP) donde cada punto representa un proceso biológico, con el tamaño correspondiente al recuento de genes asociados y el color reflejando el valor p ajustado de significancia. (B) GSEA Plot que muestra la puntuación de enriquecimiento acumulativa para la vía de regulación de la contracción del músculo liso a lo largo de la lista de genes clasificados por su expresión diferencial. (C) Dotplot de KEGG Pathways que ilustra vías metabólicas y de señalización alteradas, con el tamaño del punto indicativo del número de genes implicados y la coloración basada en la significancia estadística del enriquecimiento.
Figura 17: Análisis de enriquecimiento de vías biológicas basado en datos de expresión génica de hembras tratadas versus control. (A) Dotplot de Gene Ontology Biological Processes (GO:BP) donde cada punto representa un proceso biológico, con el tamaño correspondiente al recuento de genes asociados y el color reflejando el valor p ajustado de significancia. (B) GSEA Plot que muestra la puntuación de enriquecimiento acumulativa para la vía de regulación de la contracción del músculo liso a lo largo de la lista de genes clasificados por su expresión diferencial. (C) Dotplot de KEGG Pathways que ilustra vías metabólicas y de señalización alteradas, con el tamaño del punto indicativo del número de genes implicados y la coloración basada en la significancia estadística del enriquecimiento.

Es crucial destacar que, en las hembras, se observó la desregulación de 254 genes más que en los machos, lo que refleja una mayor sensibilidad a la sepsis inducida por LPS (Figura 18). Este fenómeno subraya la importancia de considerar las diferencias de sexo en la respuesta biológica a tratamientos y condiciones patológicas. La identificación detallada de estos genes proporciona una base sólida para investigar los mecanismos subyacentes a estas diferencias de respuesta y podría facilitar el desarrollo de estrategias terapéuticas más personalizadas y efectivas, tomando en cuenta las particularidades biológicas de cada sexo biológico.

Figura 18: Barplot que compara el número de genes significativos identificados en machos y hembras en respuesta al tratamiento, así como la diferencia neta entre estos dos grupos. Las barras azul y naranja representan la cantidad de genes significativamente alterados en machos y hembras, respectivamente, mientras que la barra verde indica la diferencia absoluta en el número de genes desregulados entre los sexos.
Figura 18: Barplot que compara el número de genes significativos identificados en machos y hembras en respuesta al tratamiento, así como la diferencia neta entre estos dos grupos. Las barras azul y naranja representan la cantidad de genes significativamente alterados en machos y hembras, respectivamente, mientras que la barra verde indica la diferencia absoluta en el número de genes desregulados entre los sexos.

Comparación machos tratados y hembras tratadas

Mediante el análisis de volcano plot y heatmap, hemos identificado una serie de genes que exhiben una expresión diferencial significativa entre machos y hembras (Figura 19A y B).

Figura 19: Análisis comparativo de expresión génica entre tratamientos en machos y hembras. (A) Volcanoplot que visualiza la diferencia en la expresión génica, con puntos de color que representan el estado de regulación: morado para genes subexpresados, amarillo para sobreexpresados y gris para aquellos no significativos, basados en la magnitud del cambio de expresión (log2 Fold Change) y la significancia estadística (-log10 p-valor). (B) Heatmap que contrasta los patrones de expresión génica entre machos y hembras tratados, utilizando una paleta de colores que varía desde el púrpura para la baja expresión hasta el verde para la alta expresión, con las muestras organizadas por similitud en el perfil de expresión.
Figura 19: Análisis comparativo de expresión génica entre tratamientos en machos y hembras. (A) Volcanoplot que visualiza la diferencia en la expresión génica, con puntos de color que representan el estado de regulación: morado para genes subexpresados, amarillo para sobreexpresados y gris para aquellos no significativos, basados en la magnitud del cambio de expresión (log2 Fold Change) y la significancia estadística (-log10 p-valor). (B) Heatmap que contrasta los patrones de expresión génica entre machos y hembras tratados, utilizando una paleta de colores que varía desde el púrpura para la baja expresión hasta el verde para la alta expresión, con las muestras organizadas por similitud en el perfil de expresión.

Entre los genes más destacados se encuentran Lima1 (LIM Domain And Actin Binding 1), que está implicado en la reorganización del citoesqueleto; Kif4 (Kinesin Family Member 4A), que juega un papel en el transporte intracelular y la segregación de cromosomas durante la división celular; Trpc4ap (Transient Receptor Potential Cation Channel Subfamily C Member 4 Associated Protein), relacionado con la señalización de calcio; Kpna1 (Karyopherin Subunit Alpha 1), que está involucrado en el transporte nuclear; Ano7 (Anoctamin 7), implicado en procesos de señalización celular y Cdc6 (Cell Division Cycle 6), esencial para la iniciación de la replicación del ADN (Figura 20).

Figura 20: Count plot de los conteos de expresión para los seis genes principales diferencialmente expresados entre tratamientos en machos y hembras, Lima1, Kif4, Trpc4ap, Kpna1, Ano7 y Cdc6. Los puntos morados indican las hembras tratadas y los amarillos los machos tratados. La posición vertical de cada punto refleja el conteo normalizado de expresión de cada gen, permitiendo la comparación entre las condiciones de tratamiento por sexo.
Figura 20: Count plot de los conteos de expresión para los seis genes principales diferencialmente expresados entre tratamientos en machos y hembras, Lima1, Kif4, Trpc4ap, Kpna1, Ano7 y Cdc6. Los puntos morados indican las hembras tratadas y los amarillos los machos tratados. La posición vertical de cada punto refleja el conteo normalizado de expresión de cada gen, permitiendo la comparación entre las condiciones de tratamiento por sexo.

Un análisis complementario mediante un diagrama de Venn permite cuantificar la distribución de los genes diferencialmente expresados entre machos y hembras enfermos. Se identificaron 1134 genes exclusivos en machos, lo que representa el 23.3% del total de genes diferencialmente expresados, mientras que en las hembras se encontraron 1391 genes únicos, constituyendo el 28.5% del total. La intersección entre ambos sexos comprende 2351 genes, representando el 48.2% del total, evidenciando una respuesta común considerable al tratamiento (Figura 21). Este dato resalta la existencia de una base común en la respuesta al LPS, aunque también se observan diferencias significativas que reflejan respuestas específicas de cada sexo.

Figura 21: Diagrama de Venn que muestra la sobreposición y la distribución exclusiva de genes diferencialmente expresados en machos y hembras. El círculo morado representa genes exclusivos de machos, el círculo turquesa representa genes exclusivos de hembras, y la intersección indica genes compartidos por ambos sexos. Los porcentajes expresan la proporción de genes únicos y compartidos en relación con el total de genes diferencialmente expresados.
Figura 21: Diagrama de Venn que muestra la sobreposición y la distribución exclusiva de genes diferencialmente expresados en machos y hembras. El círculo morado representa genes exclusivos de machos, el círculo turquesa representa genes exclusivos de hembras, y la intersección indica genes compartidos por ambos sexos. Los porcentajes expresan la proporción de genes únicos y compartidos en relación con el total de genes diferencialmente expresados.

Este patrón de expresión refleja que, si bien la principal influencia en la expresión génica parece ser el tratamiento con LPS, las diferencias entre machos y hembras también son notables, sugiriendo que la variabilidad genética asociada al sexo contribuye a modulaciones específicas en la respuesta molecular. Esto es particularmente relevante en el contexto de una medicina personalizada, donde comprender estas diferencias puede permitir la optimización de tratamientos basados en características genéticas individuales, incluyendo el sexo del paciente.

Single Cell RNA-Seq

En nuestro estudio, el análisis scRNAseq se centró exclusivamente en modelos murinos machos, a fin de desentrañar las respuestas celulares específicas involucradas en la sepsis. Para ello, se aislaron células positivas para CD11b (indicativas de macrófagos y microglia) y CD45 (marcador de células linfoides), utilizando estos marcadores para diferenciar el estado de activación y el tipo celular dentro del contexto inflamatorio de la sepsis. En este caso, las muestras a comparar corresponden a microglía extraída de ratones control y ratones tratados con LPS, todos machos.

La selección de células positivas tanto para CD11b como para CD45 nos permitió enfocarnos en macrófagos de microglia activados, células que desempeñan roles cruciales en la respuesta inmunitaria y en la mediación de la neuroinflamación dentro del SNC. Estos macrófagos activados son indicativos de una respuesta inmune intensificada y juegan un papel fundamental en la patogénesis de la sepsis, facilitando tanto la defensa contra patógenos como la potencial contribución al daño tisular asociado.

Por otro lado, las células que presentaban positividad exclusiva para CD11b fueron consideradas como microglia en un estado de activación menos pronunciado. Estas células, aunque participan en la respuesta inmune, se encuentran en un estado de activación menos agresivo comparado con los macrófagos microgliales activados, lo que sugiere una respuesta inmunitaria más modulada o en etapas iniciales de activación.

Finalmente, las células exclusivamente positivas para CD45 se interpretaron como componentes restantes del compartimento linfoides, englobando una variedad de células inmunitarias implicadas en la respuesta a la sepsis. Esta fracción representa una gama diversa de células inmunológicas, incluyendo linfocitos T, B, y otras células linfoides, que participan en la coordinación y ejecución de la respuesta inmunitaria adaptativa y que pueden jugar un papel en la progresión de la neuroinflamación y neurodegeneración en el contexto de la sepsis.

Metodología

Preparación de Datos y Entorno de Análisis

Antes de iniciar el análisis, se requiere la configuración adecuada del entorno de trabajo y la preparación de los datos iniciales. El análisis comienza con archivos FASTQ procesados por la plataforma de bioinformática del IMIB utilizando el software Cell Ranger, un método estándar en la industria para el procesamiento de datos de secuenciación de células únicas obtenidos mediante tecnologías como 10x Genomics. Los archivos resultantes, específicamente barcodes.tsv, features.tsv y matrix.mtx.gz, contienen la información esencial para la creación de objetos de análisis y son el punto de partida de nuestro estudio.

Configuración del Entorno en R

La primera etapa del análisis en R incluye la configuración del entorno y la carga de paquetes necesarios. Se utilizan diversas bibliotecas de R, que facilitan desde la manipulación de datos hasta visualizaciones avanzadas. Entre las más destacadas están: - Seurat: esencial para análisis de secuenciación de células únicas. - dplyr, ggplot2, tidyverse: para manipulación y visualización de datos. - patchwork, cowplot: para la combinación y organización de gráficos. - sctransform, ReactomeGSA, topGO: para normalización y análisis de enriquecimiento. - SingleCellExperiment, SingleR, celldex: para manejo de experimentos de células únicas y anotación de tipos celulares.

Carga y Preparación de Datos

Utilizando la función Read10X de Seurat, se cargan los datos desde los directorios especificados. Para cada conjunto de datos (Control y LPS), se crean objetos Seurat a partir de las matrices de expresión. Estos objetos son fundamentales para el manejo eficiente de los datos y facilitan la realización de análisis subsiguientes.

Control de Calidad

El análisis comienza con la evaluación de la calidad de las células, utilizando métricas como el número de genes detectados y el porcentaje de ADN mitocondrial. Se aplican filtros para eliminar células de baja calidad que puedan afectar los resultados del análisis.

Normalización y Identificación de Características Variables

Los datos son normalizados para eliminar variaciones técnicas. Posteriormente, se identifican características (genes) que muestran alta variabilidad, lo cual es crucial para las etapas de reducción de dimensionalidad y agrupación.

Reducción de Dimensionalidad y Agrupación

Se utilizan técnicas como PCA (Análisis de Componentes Principales) para reducir la dimensionalidad de los datos. Este paso es seguido por métodos de agrupación para identificar grupos o clústeres de células con perfiles de expresión similares.

Integración

El análisis culmina con la integración de todos los datos y resultados. Se realizan visualizaciones avanzadas como UMAP y t-SNE para representar las relaciones entre las células y clústeres. Además, se utilizan gráficos como heatmaps y dotplots para representar los datos de una manera que facilite la interpretación y presentación de los resultados.

Análisis de Marcadores Celulares

Para cada clúster identificado, se buscan genes que se expresan de manera única, lo que ayuda a caracterizar los tipos celulares dentro de los clústeres.

Análisis Comparativo

Se comparan las células bajo diferentes condiciones experimentales (Control vs. LPS) para identificar cambios significativos en la expresión génica que puedan ser indicativos de respuestas biológicas al tratamiento.

Análisis de Enriquecimiento

Para los genes de interés, se realizan análisis de enriquecimiento utilizando herramientas como ReactomeGSA y topGO para identificar rutas biológicas o procesos celulares relevantes.

Análisis

Los resultados obtenidos en el análisis de single-cell RNA-seq revelan aspectos cruciales sobre el impacto biológico del tratamiento con LPS en ratones, principalmente evidenciando una activación significativa del sistema inmunitario y cambios en la composición celular que podrían ser indicativos de una respuesta inmunitaria compleja o incluso de procesos patológicos similares a los observados en enfermedades neurodegenerativas.

Limitaciones en el Aislamiento Celular

Una de las principales limitaciones observadas fue el bajo número de células aisladas para el experimento, lo que podría afectar la representatividad de las poblaciones celulares y la fiabilidad de los resultados en términos de expresión diferencial de genes (Figura 22). Esta limitación es crucial, pues podría conducir a una sobrestimación de ciertos efectos o a la no detección de poblaciones celulares relevantes.

Figura 22: Violin plots que representan la distribución de las métricas de calidad de las células aisladas para control y tratamiento con LPS. nFeature_RNA indica el número de genes detectados por célula. nCount_RNA muestra el número total de transcritos por célula. percent_mito refleja el porcentaje de transcritos mitocondriales, un indicador de la salud celular.
Figura 22: Violin plots que representan la distribución de las métricas de calidad de las células aisladas para control y tratamiento con LPS. nFeature_RNA indica el número de genes detectados por célula. nCount_RNA muestra el número total de transcritos por célula. percent_mito refleja el porcentaje de transcritos mitocondriales, un indicador de la salud celular.

Caracterización de Clústeres Celulares

El análisis permitió identificar siete clústeres celulares principales, incluyendo oligodendrocitos, astrocitos, células endoteliales, monocitos, macrófagos activados, neuronas, y otras células del sistema inmunitario como células NK, células B y células T (Figura 23A). La presencia aumentada de monocitos y células del sistema inmunitario en el ratón tratado con LPS sugiere una movilización de estas células en respuesta a la sepsis. Asimismo, el incremento de oligodendrocitos podría estar vinculado a la sobreexpresión de rutas relacionadas con la diferenciación de células gliales observada en análisis anteriores (Figura 23B).

Figura 23: Análisis de poblaciones celulares mediante técnicas de aprendizaje automático y representación gráfica. (B) UMAP (Uniform Manifold Approximation and Projection) que muestra la distribución de las células en un espacio bidimensional, con el panel de la izquierda presentando todos los datos coloreados por clústeres y el de la derecha agrupado por tipo de muestra, permitiendo la visualización de la heterogeneidad celular y la identificación de subgrupos. (C) Barplot de la proporción de células por clúster, diferenciando entre control (rojo) y tratamiento LPS (azul), facilitando la comparación de la composición celular entre las condiciones experimentales.
Figura 23: Análisis de poblaciones celulares mediante técnicas de aprendizaje automático y representación gráfica. (B) UMAP (Uniform Manifold Approximation and Projection) que muestra la distribución de las células en un espacio bidimensional, con el panel de la izquierda presentando todos los datos coloreados por clústeres y el de la derecha agrupado por tipo de muestra, permitiendo la visualización de la heterogeneidad celular y la identificación de subgrupos. (C) Barplot de la proporción de células por clúster, diferenciando entre control (rojo) y tratamiento LPS (azul), facilitando la comparación de la composición celular entre las condiciones experimentales.

Análisis de Expresión Diferencial

La sobreexpresión de genes en células inmunitarias y microglía, particularmente en condiciones de estimulación con LPS, destaca genes como Atp5c1 o H2afz (Figura 24). Curiosamente, todos los genes diferencialmente expresados mostraron sobreexpresión, lo cual refuerza la idea de una activación aguda del sistema inmunitario y de respuestas celulares en el contexto de inflamación y posible daño tisular.

Figura 24: Violin plots que ilustran la distribución del nivel de expresión de genes seleccionados entre condiciones control y tratamiento con LPS. Cada plot compara la expresión en condiciones normales (rojo) con la respuesta al estímulo LPS (turquesa), con puntos individuales que representan mediciones de expresión única en células individuales.
Figura 24: Violin plots que ilustran la distribución del nivel de expresión de genes seleccionados entre condiciones control y tratamiento con LPS. Cada plot compara la expresión en condiciones normales (rojo) con la respuesta al estímulo LPS (turquesa), con puntos individuales que representan mediciones de expresión única en células individuales.

Rutas Diferencialmente Expresadas

Las rutas identificadas a través de bases de datos como Reactome, Gene Ontology (GO) y KEGG reflejan una amplia gama de procesos biológicos afectados:

  • Reactome: destacó rutas involucradas en el ciclo celular, traducción y metabolismo celular, resaltando una intensa actividad biológica en respuesta al LPS (Figura 25A).
  • Gene Ontology (GO): reveló la activación de procesos de metabolismo primarios, de proteínas y de compuestos nitrogenados, sugiriendo un estado elevado de anabolismo y catabolismo necesario para responder a la estimulación inmunitaria (Figura 25B).
  • KEGG: los genes se relacionaron con vías implicadas en la respuesta inmunitaria y procesos patológicos observados en enfermedades neurodegenerativas como Alzheimer, Huntington y Parkinson. Esto podría indicar que la activación inmunitaria por LPS comparte mecanismos moleculares con estas enfermedades neurodegenerativas.
Figura 25: Dotplots que muestran el enriquecimiento de vías biológicas basado en la secuenciación de ARN de célula única. (A) Representación de vías de REACTOME, con el tamaño de los puntos reflejando la cantidad de genes en la vía y la coloración indicando el valor p ajustado de significancia. (B) Análisis de Gene Ontology Biological Processes (GO:BP), donde la escala de puntos corresponde al recuento de genes y la intensidad de color representa el valor p ajustado.
Figura 25: Dotplots que muestran el enriquecimiento de vías biológicas basado en la secuenciación de ARN de célula única. (A) Representación de vías de REACTOME, con el tamaño de los puntos reflejando la cantidad de genes en la vía y la coloración indicando el valor p ajustado de significancia. (B) Análisis de Gene Ontology Biological Processes (GO:BP), donde la escala de puntos corresponde al recuento de genes y la intensidad de color representa el valor p ajustado.

Bibliografía

  • Catozzi C, Cuscó A, Lecchi C, Carlo ED, Vecchio D, Martucciello A, D’Angelo L, Francino O, Bonastre AS, Ceciliani F (2019) Impact of intramammary inoculation of inactivated Lactobacillus rhamnosus and antibiotics on the milk microbiota of water buffalo with subclinical mastitis. PLoS ONE 14(1): e0210204. <https://doi: 10.1371/journal.pone.0210204>

  • Rodríguez-Ramírez, K. T., Norte-Muñoz, M., Lucas-Ruiz, F., Gallego-Ortega, A., Calzaferri, F., García-Bernal, D., Martínez, C. M., Galindo-Romero, C., de Los Ríos, C., Vidal-Sanz, M., & Agudo-Barriuso, M. (2024). Retinal response to systemic inflammation differs between sexes and neurons. Frontiers in immunology, 15, 1340013. <https://doi:10.3389/fimmu.2024.1340013>

  • Salipante SJ, Kawashima T, Rosenthal C, Hoogestraat DR, Cummings LA, Sengupta DJ, Harkins TT, Cookson BT, Hoffman NG. Performance comparison of Illumina and ion torrent next-generation sequencing platforms for 16S rRNA-based bacterial community profiling. Appl Environ Microbiol. 2014 Dec;80(24):7583-91. doi: 10.1128/AEM.02206-14. Epub 2014 Sep 26. Erratum in: Appl Environ Microbiol. 2016 Jul 15;82(14):4453. PMID: 25261520; PMCID: PMC4249215.

Anexo

Figura suplementaria 1: Captura de pantalla de la salida de alineamiento de HISAT2 para la muestra 1MTreated, mostrando secuencias de ARN mapeadas contra el genoma de referencia de ratón (mm9). Cada fila representa un alineamiento individual, con columnas que indican información detallada como la posición en el genoma, cadena de alineamiento, y la secuencia del read.
Figura suplementaria 1: Captura de pantalla de la salida de alineamiento de HISAT2 para la muestra 1MTreated, mostrando secuencias de ARN mapeadas contra el genoma de referencia de ratón (mm9). Cada fila representa un alineamiento individual, con columnas que indican información detallada como la posición en el genoma, cadena de alineamiento, y la secuencia del read.
Figura suplementaria 2: Salida de Samtools stats que muestra estadísticas detalladas de la salida de alineamiento de HISAT2 para la muestra 1MTreated contra el genoma de ratón. Estas estadísticas incluyen la proporción de lecturas correctamente alineadas, la cobertura del genoma y la distribución de las lecturas en regiones genómicas específicas.
Figura suplementaria 2: Salida de Samtools stats que muestra estadísticas detalladas de la salida de alineamiento de HISAT2 para la muestra 1MTreated contra el genoma de ratón. Estas estadísticas incluyen la proporción de lecturas correctamente alineadas, la cobertura del genoma y la distribución de las lecturas en regiones genómicas específicas.